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Self-Organized Percolation Power Laws with and without 
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Abstract. Classically, percolation critical exponents are linked to the power 
laws that characterize percolation cluster fractal properties. It is found here 
that the gradient percolation power laws are conserved even for extreme gra- 
dient values for which the frontier of the infinite cluster is no more fractal. 

In particular the exponent 7/4 which was recently demonstrated to be 
the exact value for the dimension of the so-called " hull" or external perimeter 
of the incipient percolation cluster, controls the width and length of gradient 
percolation frontiers whatever the gradient magnitude. 

This behavior is extended to previous model studies of etching by a finite 
volume of etching solution in contact with a disordered solid. In such a model, 
the dynamics stop spontaneously on an equilibrium self-similar surface similar 
to the fractal frontier of gradient percolation. It is shown that the power 
laws describing the system geometry involves also the fractal dimension of the 
percolation hull, whatever the value of the dynamically generated gradient, i.e. 
even for a non-fractal frontier. 

The comparison between numerical results and the exact results that can 
be obtained analytically for extreme values of the gradient suggests that there 
exist a unique power law valid from the smallest possible scale up to infinity. 
These results suggest the possible existence of an underlying conservation law, 
relating the length and the statistical width of percolation gradient frontiers. 



1. Introduction: Gradient percolation built by diffusion or by 

ETCHING 

Spreading of objects in space with a gradient of occupation probability is most 
common. From chemical composition gradients to the distribution of plants which 
depend on their solar exposure, probability gradients exist in many inhomogeneous 
systems. In fact inhomogeneity is a rule in nature whereas most of the systems that 
physicists are studying are homogeneous as they are thought to be more simple to 
understand. In particular phase transitions or critical phenomena are studied in 
that framework, the simplest being percolation transition [J. In this work, we study 
a different situation, that of an inhomogeneous system, the gradient percolation 
situation. 
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Figure 1. The Gradient Percolation Model. Each row has an 
occupation probability p{x) ranging from 1 at the bottom to on 
the top. A site with coordinates (x, y) is occupied with probability 
p(x). The occupied and empty sites are represented respectively 
in grey and white. Apart from isolated islands and lakes, grey 
and white sites form two distinct connected regions. The marked 
sites (with thick borders) are the connected front of the occupied 
percolating cluster. 



First we study a paradigmatic situation of such an inhomogeneous system, the 
Gradient Percolation GP model [21 15] . We show how the scaling behavior, which 
usually characterizes criticality, extends down to minimal possible scales, i.e. the 
extreme gradient regime. In this regime, exact analytical analysis of the model is 
possible and given here. 

Moreover this analysis can be extended to other models, as the etching gradient 
percolation dynamical model EGP, which has recently been shown to belong to 
the same universality class as GP. In this case, the extreme gradient analysis allows 
to correctly recognize such a universality (and its limitations) without the explicit 
knowledge of the underlying dynamically generated gradient. 

Gradient Percolation GP. 

GP Definition of Gradient Percolation. The gradient percolation (GP) situa- 
tion was first introduced in the study of diffusion fronts which exists for instance as 
the consequence of a soldering process |2] |5|. It is defined in Fig.^ an d examples 
are shown in Fig. [2] The figures give examples of a random distribution of points 
on a lattice with a linear gradient of concentration in the vertical direction. 

It is a 2D square lattice of size L g x L, where each point (x, y) is occupied 
with probability p(x) = 1 — x/L g (x being the vertical direction in the figure). In 
gradient percolation there is always an infinite cluster of occupied sites as there is 
a region where p is larger than the standard percolation SP threshold p c . There is 
also an infinite cluster of empty sites as there is a region where p is smaller than p c . 
The object of interest is the GP front, the external limit (or frontier) of the infinite 
occupied cluster. 
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Figure 2. Gradient percolation (GP) front. Particles are dis- 
tributed at random with probability p(x) = 1 — x/L g . The occu- 
pied sites are in black and the last line of connected occupied sites 
is the gradient percolation front shown in light gray. This situation 
corresponds also to a diffusion situation where the front is called 
the diffusion front. Top: L = L g = 500. Bottom: L = 50 and 
L g = 12. The black window has an horizontal width equal to (jf. 
It contains approximately L g points. 

Its precise definition depends on the lattice geometry: 

• Triangular Lattice. The GP front is the connected set of sites belonging 
to the occupied infinite cluster which are first nearest neighbours to sites 
of the infinite empty cluster. 

• Square Lattice. The GP front is the connected set of sites belonging to 
the occupied infinite cluster which are first or second nearest neighbours 
with sites of the infinite empty cluster, itself defined through a first and 
second neighbours connection. 

The square lattice case is shown in grey in Fig. [5] This front is a random object 
with an average position xt, a statistical width <jf and a total length Nf. In so far 
that the GP front and the SP external perimeter (often called hull) have the same 
geometry their fractal dimension was first conjectured to be exactly equal to 7/4 
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Figure 3. Sketch of the etching dynamics in a square lattice: 
sites 2,3,5 are etched at the first time-step as their resistances are 
smaller than p(0). Consequently the number of etchant particles 
in the solution decreases by 3 units. At t = 1, the new interface 
sites than can possibly be etched are 7, 8, 10 if the solution can 
etch only the first nearest neighbours. If the solution can also etch 
second nearest neighbours in a diagonal direction the whole second 
layer can possibly be etched. 

in |2]. This result was then demonstrated heuristically by Saleur and Duplantier 
|3] and very recently it was proved mathematically by Smirnov and Werner p>] . 

Etching Gradient Percolation EGP. The same type of fractal geometry 
has been recently found in a very different physical situation. It first appeared in 
experiments [B]. It was interpreted in a model in which an etching solution is in 
contact with the initially flat surface of a disordered solid and starts to corrode its 
weakest regions. The solid surface then gets "harder" but at the same time new 
regions are discovered which contain weak elements. Often the corrosive power 
of the solution is proportional to an etchant concentration and if the etchant is 
consumed in the reaction, the corrosive power of a finite volume of solution decreases 
during the time evolution of the process. As the solid surface gets harder and 
harder, and the corroding power gets weaker and weaker, the corrosion process 
stops spontaneously in a finite time interval. 

The EGP model definition. We first recall the two-dimensional etching model 
introduced in 7 . Its schematic is shown in Fig. [3J 

• A 2D random solid is represented as a site lattice (triangular or square), 
of linear width L and, eventually, infinite depth. 

• A random number rj <E [0, 1] (extracted from the flat probability density 
function 7To(r) = 1 for r £ [0, 1]) is assigned to each solid site i, represent- 
ing its resistance to the etching by the solution, rj does not depend on 
time (quenched disorder), or on the site environment. 

• The etching solution has a volume V and is initially in contact with the 
solid through the bottom boundary (see Fig. OJ). It contains an initial 
number N et (0) of dissolved etchant molecules. 

Consequently, the initial concentration C(0) of etchant in the solution is given 
by: C(0) = N et (0)/V. Calling N et (t) the number of etchant molecules at time t, 
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Figure 4. Etching process for two intermediate times and final 
equilibrium situation. The solid is shown in grey, the solution in 
white, and the finite size solid clusters that are detached by the 
etching in black. The final solid surface is found to be fractal with 
dimension 7/4 up to a characteristic scale cr/. 



C{t) = N et (t)/V. At each time-step, the "etching power" of the solution (i.e. the 
average "force" exerted by the solution on a solid surface particle) is supposed to 
be proportional to C(t) : p(t) = TC(t). Hereafter the assumption T = 1 is made, 
without loss of generality. It implies C(t) = p(t). At time-step t, all the interface 
sites with rj < pit) are dissolved and a particle of etchant is consumed for each of 
these corroded solid sites. 

The etching evolution is depicted in Fig.|I] A brief review of the known results, 
will be given in the next section. 

The structure of this paper is the following: in the next section we recall 
definitions and known results for fractal "weak gradient" cases. Then we consider 
the extreme gradient case for GP. In that case a few exact results are given. They 
are discussed in the spirit of deciding whether or not there exists power laws for 
GP that are valid from the smallest scale to infinity. 

They are then compared with numerical results obtained for small size systems. 
Next the EGP is studied numerically for large gradients. The results are then 
extrapolated to large (fractal) systems and a remarkable agreement is found between 
extrapolated values and numerical values. Gradient percolation results are then 
compared with etching gradient percolation. In the conclusion we discuss to what 
extent the numerical results give evidence that there exists a new conservation law 
in diffusion or gradient percolation, this law being valid from the smallest scale to 
infinity. 
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2. GP and EGP Known Results 

2.1. Gradient Percolation (GP): Known Results. The early GP studies 
were focussed at finding its relations with standard percolation. Let us first recall 
the definitions. For < x < L g , nf(x) is the mean number of points of the front 
lying on the line x per unit horizontal length. It measures the front density at 
distance x. The length Nf, the position Xf and the width <jf of the front are then 
defined in terms of Tif{x) by 

, > T,x=0 Xn f( X ) 



N f = ^L^W' x f-^IT g 7T 

x=0 X-ix=O n f\ X ) 



and rf = D^(*-*/)V(*) 



f T,xLo n f( x ) 

It was found that the mean front was located at a distance where the density of 
occupation was very close to p c or p{xj) ~ p c - This was verified numerically with 
such precision that the gradient percolation method is now often used to compute 
percolation thresholds j3 El 1101 111! I12| . It was also found that: 

(1) The width <j/ depends on L g through a power law <jf oc {L g ) v ^ 1+ ^ where 
v = 4/3 is the correlation length exponent [l] in dimension d = 2 so that 
at = (L g ) 4 / 7 . The width Of was also shown to be a percolation correlation 
length. 

(2) Secondly it was found that the front was fractal with a dimension Df, 
numerically determined, close to 1.75. The front length followed a power 
law N f cx (L g ) aN with a N = (D f - l)v/(l + v). 

(3) But also, it was numerically observed that the sum of these two exponents 
was very close to 1 . If true, this meant that v/{l+v)-\-{Df—l)v/(l+v) = 1 
or Df = 1 + 1/v = 7/4. This is how it was conjectured in [2] that 
D f = 7/4. 

In that sense the ordinary GP power laws were thought to be linked to the SP 
exponent v and to the fractality of the percolation cluster hull. Up to now, these 
facts were considered to be strictly valid only in the large system limit. 

However, if true, and we know now that 7/4 is the exact value, there follows 

Dt 

an intriguing relation, namely a , is exactly proportional to L g . This means that 
the number of surface particles within the correlation length is exactly proportional 
to L g . This is particularly striking for diffusion fronts. Diffusion of particles from 
a source results in a concentration gradient and an associated GP situation. In 
that frame, the above result means that, if L g particles have diffused on a vertical 
row, there is on average the same (or a constant fraction of) number of particles 
on the correlated surface (surface content of a box with a lateral size equal to the 
statistical width). This fact seems a priori to have nothing to do with scaling, 
percolation and the thermodynamic limit. From this point of view it is possibly 
the consequence of a conservation law and if such a conservation exists, it should 
apply also for extreme gradients corresponding to L g of a few units. 

In particular it should apply to the very extreme L g = 1, 2 and 3 for which 
exact values of Xf, Nf and 07 can be calculated analytically. The corresponding 
results are the subject of section 3. 
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Figure 5. Box-counting determination of the fractal dimension 
Df of the corrosion front (e is the linear size of the box, N(e) the 
number of boxes containing at least one point of the front). The 
value of Df = 1.753 ± 0.005 is found fit for values of e ranging 
from a few lattice distances to the front width aj (in this case 
a f w 3000, i.e. 01(07) « 8). 



2.2. Etching Gradient Percolation (EGP): Known Results. Let us call 
n(t) the number of dissolved solid sites at time-step t. One can express several 
quantities through n(t), or its time-integral N(t), that is the total number of cor- 
roded solid sites up to time t. The number of etchant particles in the liquid will 
decrease as: 



(1) 



N et (t + 1) = N et (t) - n(t) = N et (0) - N(t) , 



and consequently the etching power of the solution is: 

n(t) N(t) 



(2) 



p(t + l)=p(t) 



V 



= p(0) 



V 



Note that, as p(t+ 1) < p(t), a site having resisted to etching at a certain time-step 
will resist forever. Consequently, the part of the solid surface which can be etched at 
time-step t+1 is restricted to the sites which have been just uncovered by the etching 
process at time t. We call this subset of surface the "active" part of the surface. 
After a given time-step, all the solid sites which have been previously explored by 
the solution are "passive" forever. However it may happen that "passive" sites are 
disconnected from the bulk at a later time-step if they are connected to the solid 
by weak sites. 

The model reproduces qualitatively the same phenomenology observed experi- 
mentally [B]. The dynamical evolution can be divided into two different regimes: 
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Figure 6. Behavior of af and of the average number of sites per 
column Nf /L as a function of V/ L for several sample sizes (several 
values of L). 



(1) In the first (smooth) regime, the corrosion is well directed while the front 
becomes progressively rougher and rougher. In our model this regime does 
not depend on the details of the discretization, or on the fundamental 
geometrical features of the lattice, like the embedding space dimension or 
the lattice coordination number. 

(2) In the second regime, the correlations revealed by the hardening process 
become important: the dynamics becomes locally isotropic, generating a 
fractal front. This corresponds to a critical regime, directly related to the 
percolation transition on the same lattice. 

In the case of etching the same results as for GP were found on the final front 
geometry. The fractal dimension Df of the etching front was measured (up to the 
scale <jf) using the box-counting |13| algorithm. In this way Df = 1.753 ±0.005 is 
measured (see Fig.EJ very close to the value 7/4 of GP. 

In fact it was shown |14| that in EGP the length V/L, ratio of the solution 
volume by the sample width, plays the role of L g . The width and length of the 
front were found to follow the same power laws as in GP, as shown in Fig. [5] In 
particular <jf scales as (V/L) a ", with a a 0.57, while Nf/L scales as (V/L) aN , 
with an observed value for a at ~ 0.45. 

The scaling relation a a + apj = 1 seems also fairly obeyed, but a direct in- 
spection give an even better result, as shown in Fig. [7| where afNf/L is plotted 
as a function of V/L. A power law fit gives an exponent very close to one, and 
correspondingly, a linear fit gives equivalently good results. 

As the same relations are verified by EGP, there should also exist a conservation 
law for that case. This law would stipulate that the amount of correlated final front 
is proportional to the ratio of the solution volume divided by the lateral size of the 
sample, in other words to the depth of the solution. Formulated in this manner 
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Figure 7. Behavior of <jfNf/L as a function of V/L. Data are 
fitted with a power law, resulting in an exponent very close to one, 
and with a linear law, resulting in comparable results. 

there is no reason why this property should not be verified for very small solution 
volumes, as shown in the next section. 



3. Extreme Gradients: fractal exponents without fractal 

geometry 

In this section we will show how, even in the case of extreme gradient (hence 
in absence of any fractal structure) one recovers the above power laws. 

3.1. Exact calculations for extreme GP. We present a detailed descrip- 
tion of some exact results obtained for L g = 2, 3, hence the two highest non trivial 
value for the gradients. 

3.1.1. Triangular lattice. 

Exact computations for L g = 2. We find occupied sites only on the lines x — 
and x = 1. For 1 < i < L, let Xi be the random variable which has value 1 if the \ th 
site of the line x = 1 is occupied, and otherwise. Then, the X^s are independent, 
and each takes value 1 with probability 1/2 (see Fig. |SJ. 

• Since all the occupied sites on the line x — 1 belong to the front, we have 
Nf(l) — Yli=i Xi- The mean number of points of the front on this line is 
thus E(iV/(l)) = L/2. 

• On the line x = all the sites are occupied, but a site will belong to 
the front if and only if at least one of its two neighbours on the line 
x = 1 is not occupied (see Fig. |SJ). Thus, Nf(Q) = L — J2f =1 

and the mean number of points of the front on the line x = is then 
E(iV/(0)) = L-L/4= 3L/4. 

Finally, we get 

3 1 
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L 

Figure 8. Gradient percolation on the triangular lattice with 
L g = 2 (the sites on the line x = 1 are occupied with proba- 
bility 1/2). Thick disks represent occupied sites, and the thick line 
joins the sites belonging to the GP front. 




Figure 9. Gradient percolation on the triangular lattice with 
L g = 3 (the sites on the line x — 1 are occupied with proba- 
bility 2/3, and the ones on the line x — 2 with probability 1/3). 
Thick disks represent occupied sites, and the thick line joins the 
sites belonging to the front. 



and thus we can compute the mean length Nj, the mean position Xf and the width 
Of of the front, and find: 

^1 = - = 1.25 , at = ^ = 0.4898979... , 

L 4 1 5 

x f = - = 0.4 and p(x f ) = 1 - ^ = 0.8 . 

5 J Lg 

Exact computations for L g — 3. For 1 < i < L, let be the random variable 
which has value 1 if the i th site of the line x = 2 is occupied, and otherwise. 
Then, the X^'s are independent, and each takes value 1 with probability 1/3. In 
the same way, we also define Yi as the random variable which has value 1 if the i th 
site of the line x = 1 is occupied, and otherwise. The Y^s are also independent, 
and each takes value 1 with probability 2/3 (see Fig. EJ. 

• A site on the line x = 2 will belong to the front if it is occupied and if 
it does not belong to an isolated cluster (i.e. a cluster of occupied sites 
on the line x = 2, surrounded by empty sites, see Fig. EJ) . We can write 
this in terms of X$ and Yi in the following way: Nf(2) = X\ + ... + Xl — 

EtlEnM^-^-l^-^+n-lfl-^+n+lJtl-y.l-tl-y.+n). The 

mean number of points of the front lying on the line x = 2 is thus 

E(Nf(2)) = -L-Lx-x-xy^n—= —L. 

3 9 3-^9" 16 

TI>1 
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• A site on the line x = 1 will belong to the front if it is occupied and if: 
either at least one of its two neighbours on the line x — 2 is empty, either 
these two neighbours are occupied but the site is connected to the "sea" 
on its left or on its right (see for example, on Fig. |5J the third occupied 
site from the left on the line x = 1). Thus, 

L L 

N f (l) = YiQ- - XiXi-i) + Y i x i x i-i(Si + S[ - S^), 

i=l i=l 

where Si = Y^ntii 1 ~ ~ K i _„)X 1 _ 2 ---^i-n(l - X-n-i) is a bi- 

nary (0 or 1) random variable which has value 1 when the site i (on the 
line x = 1) is connected on its left to the "sea". S- is the analogous of 
Si for the connection on its right. We notice that S^ and Si are indepen- 
dent, identically distributed, and their mean value is E(Si) = E(5j') = 
E„>i(l/3)"(V3)"- 1 x2/3 = l/4. 

The mean number of points of the front lying on the line x — 1 is thus 

E(N f (l)) = ~(1 - -)L+ - x -12 x - - 4)i = ~L. 

• A site on the line x = will not belong to the front if: else, its two 
neighbours on the line x = 1 are occupied, or else, one of these two 
neighbours is empty but the site is not connected to the "sea" (this is for 
example the case of three consecutive sites of the line x = on the right 
of Fig.©. Thus, we have: JV>(0) = L — Y.ti Y i Y i+i - Eti £„>i(™ + 
l)Yi(l - Y i+1 )...(l - Y i+n )Y i+n+1 X i X i+ - L ...X i+n . Then, the mean number 
of points of the front lying on the line x = is 

MN f (0)) = L L — L x - x - VVn + 1)— = — -L. 

y n " 9 9 3 ^ v ' 9™ 432 

n>l 



To summarize, we get 



223 5 5 

nf(0) — , nr(l) = - and n#(2) = — . 

n ' 432' n ' 8 n ' 16' 

and then we can compute the mean length Nt, the mean position Xf and the width 
<jf of the front, and find: 

*k = 15_ = 1.453704... ,a f = ^ = 0.7419084... , 
L 108 1 157V2 

x f = ^ = 0.8598726... and p(x f ) = 0.7133758... . 

3.1.2. Square lattice. The situation for the square lattice is a little different 
from the situation for the triangular lattice, because of the problem of the choice 
between 4-connexity and 8-connexity. In our case, the definition of the front will 
be the following: the front is the set of occupied sites which are connected through 
occupied 4-neighbours to the line x = 0, and have an empty 8-neighbours which 
belongs to the empty 8-connected component of the line x = L g . 




Exact computations for L g — 2. The computations for L g = 2 in the case of the 
square lattice are very similar to the ones in the triangular case. All the occupied 
sites of the line x = 1 belong to the front, and a site on the line x = belongs to 
the front if at least one of is three neighbours on the line x = 1 is empty. Thus, we 
get 

17 1 

= 1 - ^3 = g and n f0) = 2' 
and then we can compute the mean length Nf, the mean position Xf and the width 
07 of the front: 

^- = 11 = 1.375,^ = ^=0.4810457... 
L 8 1 11 

x f = Jl= 0.36363636.... and p(x f ) = 0.81818181 

Exact computations for L g — 3. For L g = 3, the computations are similar to 
the ones in the case of the triangular lattice, but a little more complicated since a 
site has more neighbours, and thus more different geometric configurations have to 
be considered. But it can be done, and we finally get: 

8029 . , 47 , 13 
n ' (0)= 11664' n ' (1)= 72 md n ' (2)= 48' 
and for the mean length Nf, the mean position Xf and the width 07 of the front, 
we find: 

N f 9401 9V576049 
-L = — — = 1.611968... ,a f = — — = 0.7266046.... 
L 5832 ' 1 9401 

6966 

x f = = 0.7409850... and p(x f ) = 0.7530050 

1 9401 v " 

3.2. Numerical results. The problem is then to compare the numerical GP 
laws to these exact values. As will be shown, the numerical results verify the above 
power laws with such precision that the question arises of the existence of a simple 
mathematical power law extending from L g = 1 to infinity. To try to answer 
this question we proceed in two steps. First we test these laws on the numerical 
results obtained for L g between 4 and 50 for the square and triangular lattices 
by searching the best numerical power laws followed by the width. Considering 
arbitrary exponents a between 1.6 and 1.9, we study aj as a function of L g between 
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FIGURE 11. Numerical results for the square lattice. Left: the 
circles resp. diamonds represent respectively a 7 / 4 and (Nf/L) 7 / 3 . 
Right: the stars represent Nf x <Jf/L. The lines are the ones 
obtained by linear fit on the data for L g — 4 to 50. 



square lattice 


L g = l 


Lg = 2 


L g = 3 


exact a f 





0.4810 


0.7266 


a f (4-50) 


-0.13 


0.47 


0.72 


a f (4-5) 


-0.106 


0.478 


0.726 


Sa f (4-5) 


0.025 


0.006 


0.002 


exact Nf/L 


1 


1.3750 


1.6120 


Nf/L (4-50) 


1.24 


1.48 


1.68 


N f /L (4-5) 


1.109 


1.393 


1.615 


S(N f /L) (4-5) 


0.017 


0.009 


0.004 



Table 1 . Comparison between exact and extrapolated results for 
the square lattice. The data (4 — 50) resp. (4 — 5) correspond to 
extrapolated values from the respective ranges (4 — 50) resp. (4—5) 
(see text). 5 is the confidence interval. 



4 and 50. For each a value, there is a best line aj = a a (L g + b a ) fitting the 
numerical aj. The introduction of the term b a is justified by the fact that when 
a power law is verified for large systems it includes always the possibility that a 
small (as compared to the system size) term could contribute but in a negligible 
manner. But here the size itself is small or very small. On the other hand, one 
should remark that for L g = 1 the width is strictly so that some negative value 
of b a should be present. In the next step, the mean error d(a), defined by d(a) 2 — 
(1/47) X)z° s =4( cr /(^s) Q — a a(L g + b a )) 2 , is measured numerically as a function of a. 
The results are shown in Fig. 1101 There is a clear minimum for a 1.75, showing 
that this exponent gives the best power law fit. Once the best fit with the empirical 
data is made one has the best values for the parameters a and b: a = 0.297 and 
b = —1.094. Note that b should be strictly equal to —1 in order to obtain a null 
width for the trivial case L g = 1. 
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triangular lattice 


L g = l 


Lg = 2 


L g = 3 


exact <jf 





0.4899 


0.7419 


Of (4-50) 


-0.15 


0.48 


0.74 


a s (4-5) 


-0.151 


0.480 


0.740 


6v f (4-5) 


0.019 


0.006 


0.002 


exact Nf/L 


1 


1.2500 


1.4537 


N f /L (4-50) 


1.16 


1.36 


1.46 


N f /L (4-5) 


1.014 


1.262 


1.457 


S(N f /L) (4-5) 


0.015 


0.008 


0.004 



Table 2. Comparison between exact and extrapolated results for 
the triangular lattice. The data (4 — 50) resp. (4 — 5) correspond 
to extrapolated values from the respective ranges (4 — 50) resp. 
(4 — 5) (see text). 6 is the confidence interval. 



Another verification of the extreme GP power laws can be obtained from the 
study of the front length or of the quantity (Nf/L) 7 / 3 as a function of L g . In 
Fig. 1111 (left plot), the diamonds represent the values of (Nf/L) 1 / 3 and the best 
linear fit has equation Y = c(L g + d) with c = 0.843 and d = 0.959. This shows 
indeed that the exponents 4/7 and 3/7 can be used down to the steepest gradients 
for which the frontier is no more fractal. 

We also study the quantity Nf x Of/L as a function of L g : the result is shown 
on Fig. ^2 (right plot). The best linear fit (on the values from L g = 4 to 50 has 
equation Y = c(L g + d) with c = 0.465 and d = -0.278. 

One can then extrapolate the a* values to the case L g = 1, 2 and 3. The results 
are given in Table for both lattices. One observes that the numerical extrapo- 
lations correspond to the exact values with an apparent good precision. However 
no firm conclusion can be drawn without discussion of the numerical uncertainties. 
The values shown on Fig.^Jare averaged over 100 trials on a length L = 5.10 5 . As 
the fit occurs through a power law it is difficult to give a confidence interval for the 

coefficients a and b (obtained from a least square linear regression on the values of 

7/4, 
a/ ). 

In order to obtain a better control on the numerical precision of a and b we 
made extensive computations of the two cases L g = 4 and L g = 5 with 100 trials on 
a length L — 5.10 5 . Doing so we obtain the mean values with their standard devia- 
tion: cr/(4) 7 / 4 = 0.8658±0.0009andforcr/(5) 7 / 4 = 1.1610±0.0013. Thus if we com- 
pute the equation of the line a(L g +b) which interpolates the two points (4, <r/(4) 7 / 4 ) 
and (5,cr/(5) 7 / 4 ), we obtain a = 0.2952 ± 0.0022 and b = -1.066 ± 0.041. This last 
result shows that the value —1 is compatible with b and its statistical error. Given 
the numerical values for L g = 4 and 5, we can also get extrapolated values for of 
for L g smaller, together with their confidence interval. The result (see Tables Q] 
andEJ is that the predicted values are very close to the exact ones. For (Nf/L) 7 / 3 , 
in the same way we obtain a linear interpolation of the values for L g = 4 and 5, 
with coefficients c = 0, 893 ± 0.014 and d = 0.427 ± 0.074. 

A comparison of the main fitting results for the two geometries (square and 
triangular) is shown in Table 01 
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numerical results 


square lattice 

100 x 500000 


triangular lattice 

100 x 500000 


er^ 4 fit for 
data Lg = 4 — * 50 


a = 0.297, 
b= -1.094 


a = 0.315, 
b = -1.113 


(7V//L) 7 / 3 fit for 
data L# = 4 — > 50 


c= 0.843, 
d = 0.959 


c = 0.636, 
d= 1.214 


a f x Nf/L fit for 
data Lg = 4 — > 50 


5 = 0.465, 
d= -0.278 


c = 0.426, 
d= -0.198 


exact linear fit from 
a 7 / 4 (2) and a) 1 ' 4 (3) 


a = 0.2940, 
& = -1.0548 


a = 0.3062, 
b = -1.0632 



Table 3. Numerical results of the linear fit for the square and 
triangular lattices. The data are obtained from the averaged over 
100 trials with a length L = 500000. 



3.3. Extreme Gradients for EGP. To check the extreme gradient regime 
for EGP numerical simulations have been performed for large systems (L > 5000) 
in a wide range of etchant volumes V. As discussed in 14 , the ratio L/V is 
proportional to the self-established gradient generated by the dynamical process, 
at least in the small gradient regime. Here we show how it is possible to check the 
whole range of gradients (from small to extreme) without explicitly knowing the 
exact value of the gradient. 

Let us assume that there exists such a parameter g, or equivalently L g = l/g, 

4/7 3/7 

such that af oc L g and Nf cx L g . In such a case, one should have, along the 
lines discussed above: 

(3) =a{^-) V \b 

In order to check the validity of the above relation on the entire gradient range, 
we proceeding in the following manner. First, we use relation © to determine the 
parameter A and B fitting the data on an extreme gradient (non fractal) range. 
Once the parameters A and B have been determined, we use, we use relation [3] to 
check the agreement with the numerical results obtained for small gradients, in the 
fractal regime. In Fig. El the outcome of such a procedure is shown, in the case of 
a triangular lattice. In the first graph (top left) the fitting region is shown. The 
data correspond to the simulation of more than 100 realizations of systems sized 
L > 10000 for V ranging from V = L to V w 3L. The resulting values of Nf/L 
correspond, in the GP model, to L g ranging approximatively from 1 to 70. The 
result of the fit is A = 0.58804 and B = -0.89278. The following four graphs (from 
left to right in the first and then in the second rows) display the values predicted 
by the linear relation (straight line) compared with the simulation results (for 
a system sized L = 5000). Each graph corresponds to a scale ten times larger than 
the previous graph. Finally, the last graph (bottom right) shows the whole range 
of simulation results together with the linear relation © , in a log- log scale. 

3.4. Comparison between Gradient Percolation and Etching Gradi- 
ent Percolation. The process of creation of the interface being different, the ques- 
tion arises of the comparison between the two geometries. For this we compare the 
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Fitting region ...larger... ...larger... 




Figure 12. E.G.P numerical scaling extrapolated from the ex- 
treme gradient regime. In the first graph (top-left) the data are 
fitted against a linear relation. Then the predicted values for larger 
ranges are compared with simulations. The last graph (bottom- 
right) displays the whole range of data, in a log-log plot. 



extreme gradient results of the two models. In Fig.^Ilwe show the numerical results 
for GP and EGP for both lattices. The result of the linear fit, with the relation © 
are displayed in Tabled 

First note that, as expected, the slopes are different for the two lattices. More 
interestingly, given a lattice geometry, the GP model and the EGP model have 
slightly different coefficients A and B. This shows that GP and EGP are not 
strictly equivalent, even if they belong to the same universality class (with respect 
to exponents). The differences between GP and EGP can be explained in the 
following way. In EGP the probability of invasion by the solution is applied on the 
previous irregular frontier while in GP the probability changes row by row. In that 
sense in EGP the front penetration may be slightly increased as compared to GP. 

Another quantity that displays a difference between the two models, is the 
depth (i.e. x-coordinate) distribution of the front sites, as shown in Fig. El Note 
that, at odds with the GP case, the EGP is asymmetric. To understand such an 
asymmetry one has to recall the mechanism that creates a gradient in EGP. There, 
at each time t there is an etching power p(t) and the front is at an average distance 
x(t). This means that p depends implictly on x, through the dummy variable 
t. Consequently, there is an effective "self-generated" gradient, but there is no 
reason why this gradient should be constant as it is the case for standard GP. In 
particular in the last period of etching, p varies only weakly below p c |14j while the 
penetration still increases. This means that the gradient corresponding to the last 
period is smaller, inducing a wider distribution for the front. Alternatively if GP 
was studied with a concave p(x) (decreasing gradient) the frontier would naturally 
be asymmetric. 
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(N/L) 



Figure 13. Extreme Gradients: GP (empty symbols) versus EGP 
(filled symbols) for square and triangular geometries. 



GP triangular 

Lg=4000, L=20000 <p>=0.5 

i 1 i 1 r~ 



GP square 

Lg=3000, L=20000 <p> = 0.59259 

n i 1 i 1 r 
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2000 
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000 



EGP triangular 

L=20000, V=1.4xl0 7 



1300 



EGP square 

L=20000, V=1.25xl0 7 



1400 




100 200 300 400 500 600 



700 



Figure 14. Depth distribution for the sites of the final front. Note 
the asymmetry for the EGP model, in contrast to the symmetric 
GP distributions. 



4. Discussion and Conclusion 

We have seen above that the detailed study of the extrapolation from the best 
numerical values corresponding to L g = 4 or L g = 5 was, for the square lattice, 
compatible with the exact values within numerical uncertainties. At this point 
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square lattice 


A 


B 


EGP 


0.40808 


-0.67458 


GP 


0.35142 


-0.58746 


triangular lattice 


A 


B 


EGP 


0.58804 


-0.89278 


GP 


0.49459 


-0.70487 



Table 4. Comparison between GP, EGP. The table shows the 
results of the fit performed with the relation J3J on the data shown 
in Fig. HI 



one could conclude that the numerical results are compatible with the existence 
of a single mathematical law for the width dependence, this law being valid from 
L g = 1 to infinity. Note that the quality of the random number generator could 
intervene in the purely numerical results and uncertainties, but when we tested the 
case L g = 2 or L g = 3 we did obtain numerically the exact results. Still, as long 
as a mathematical proof has not been given, it is not possible to conclude on the 
exact values of the coefficients a and b. 

The question of a unique mathematical power law can also be addressed from 
the exact results only. One can note that, remarkably, each time one tries to fit an 
observable y (as (iV//L) 7 / 3 or er 7 / 4 ) using a function of L g like: 

V = a ( L g + b ) 

the values of the term b in the fit of the width are close to —1 (see Table 0J. If 
the above law exists from L g = 1 to infinity, it suggests that the real value of b is 
exactly -1 as the width is null in the trivial case L g — 1. As we have exact values, 
one can compute the equation of the line y = a(L g + b) defined by the two points 
(2,cr/(2) 7 / 4 ) and (3, fj/(3) 7 / 4 ). One obtains a = 0.294 and b = -1.055. These 
values are close to the values obtained from the above numerical fit (a = 0.297 and 
b = —1.09), and here again b is close but not equal to —1. Why is there a small 
mismatch with the simplest law? The answer to this question is two-fold. 

(1) It is possible that the exact power law is not valid for L g = 1 or both 
for L g = 1 and 2. These cases could be "anormal" as in these cases 
the perimeter and the accessible perimeter are the same (i.e. there is no 
Grossman- Aharony effect |15| ). It is then possible that these two cases 
cannot be explained by the same mathematical law. 

(2) The small discrepancy could be related to the fact that a slightly different 
definition of the interface, leading to very similar results for large values 
of L g (where universality enters), can give slightly different exact values 
for Lg = 1,2,3. Our frontier definition considers only the occupied sites. 
It gives to these sites a privilege role whereas one should also consider the 
frontier of the empty cluster. In fact this is not new in percolation studies 
[51 where it was shown that the barycenter between the frontier of 
the occupied cluster and the frontier of the empty cluster was a more nat- 
ural object. It notably permitted better computations of the percolation 
threshold. We have studied the statistical width of the local barycenter 
which can also be computed exactly for L g — 2 or 3. The results show the 
same behavior as described above i.e. a value of b close but not equal to 
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— 1. The question remains then open to define the nature of the geomet- 
rical object which would really display a b value exactly equal to —1. The 
answer to that question would certainly have interesting consequences for 
the SP problem itself. 

In summary, it has been shown that the classical power laws of gradient per- 
colation can be extended to extreme gradients with the same fractal exponents 
although the systems present no fractal geometry. Several comments can be drawn 
on these results. 

First, from the purely mathematical point of view, our results suggest that 
there exists a conservation law which stipulates that the length of the correlated 
frontier is strictly proportional to the gradient length. This hypothesis remains to 
be proved rigorously. 

Secondly, in a wider theoretical frame, the fact that the exponents 4/7 and 3/7 
are valid down to the smallest L g values (or the steepest gradients) suggests that 
these exponents play the same type of role here that the exponent 1 /2 intervening 
in the fluctuations of the sum of independent identical random variables. In that 
last case the exponent applies to any number of random variables starting from 1, 
2 or 3 up to infinity. The exponents 4/7 and 3/7 may then play a more important 
role than critical exponents which only exist in the thermodynamic limit. 

Third, the fact that the same exponent has been found for the square and the 
triangular lattice (and in two different models) marks a universal behavior that, 
in principle, is quite unexpected here. Up to now, the main physical argument 
to explain universality, i.e. the independence of exponents from the microscopic 
details of the model, was based on the long range of the correlations which appear 
in the proximity of a critical point. In the extreme gradient case, on the contrary, 
the (geometrical) correlations of the interface are so small that one cannot justify 
the observed universality. 

Finally, our study suggests an intrinsic method to determine whether a given 
rough interface belongs to gradient percolation, without knowledge of the gradient 
and whatever the width of the interface. This can be very helpful to understand 
the properties of diffused contacts between materials. 
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